General global pattern of diversity with elevation
Bird community data
Resolution: 100 m elevational bands
Figure 1: Locations of the 46 elevational gradients used in the study.
Four diversity metrics
Quadratic regression with elevation
Classify curve shapes
Four diversity metrics
Quadratic regression with elevation
Classify curve shapes
Consistency across mountains & metrics
1. Elevational patterns of diversity
Figure 2: Patterns of bird functional and phylogenetic diversity along elevational gradients.
2. Deterministic strength
2. Deterministic strength
Figure 3: Global trends of standardised effect sizes of diversity metrics against standardised elevation, as predicted from quadratic regressions and a global quadratic mixed-effects model.
3. Latitudinal effects
3. Latitudinal effects
Figure 4: Effect of latitude on deterministic strength of elevational trends in diversity.
# iterate over each mountain in the dataset:# fit LM of each metric against linear and second-order polynomial terms of Elevation# extract intercept and coefficients for both Elevation terms# aggregate all these values over all mountainslist_mount <-unique(data_dryad$Mountain)metric_cols <-c("FRic", "FDis", "PD", "MPD")to_iter <-expand.grid(mountain = list_mount, metric = metric_cols)quad_reg <- purrr::map2(to_iter$mountain, to_iter$metric, ~ { data_mount <- data_dryad %>%filter(Mountain == .x) model_formula <-as.formula(glue("{.y} ~ elevation + I(elevation^2)")) model_mount <-lm(model_formula, data_mount)# predict only for elevations present in each mountain to_pred <- data_mount %>%distinct(Mountain, elevation)# predict model_pred <-predict(model_mount, newdata = to_pred, type ="response") pred <-bind_cols(to_pred, pred = model_pred)# store p-value of model fit p <- broom::glance(model_mount) %>%pull(p.value)# bind all outputs togethert(model_mount$coefficients) %>%as_tibble() %>% magrittr::set_colnames(c("intercept", "coef_linear", "coef_poly")) %>%bind_cols(tibble(mountain = .x, metric = .y, pval = p)) %>%left_join(pred, by =c("mountain"="Mountain"))}) %>%list_rbind() %>%relocate(metric, mountain) # bring metric and mountain cols to the start
Plotting predicted patterns
Code
quad_reg %>%ggplot(aes(x = elevation, y = pred, col = mountain)) +geom_line() +facet_wrap(~ metric, ncol =1, scales ="free_y", strip.position ="left") +labs(x ="Elevation", y ="") +# scale_col_pal + # only makes sense to colour scheme after classifying into 8 categstheme(panel.grid.major.x =element_line(colour ="grey", linetype ="dotted"),# match strip format in paperstrip.background =element_rect(colour =NA),strip.placement ="outside",strip.text =element_text(size =10),# remove colour keylegend.position ="none")
# plot frequency bar graph of different trend typesquad_reg_class %>%mutate(trend =factor(trend, levels =c("Increasing", "Decreasing", "Mid Peak", "Mid Valley","Low Peak", "Low Valley", "High Peak", "High Valley", "NS"))) %>%group_by(metric) %>%count(trend) %>%ggplot() +geom_col(aes(x = trend, y = n, fill = trend)) +facet_wrap(~ metric, ncol =1, strip.position ="left") +labs(x ="Trend type", y ="") + scale_fill_pal +theme(legend.position ="none",# match strip format in paperstrip.background =element_rect(colour =NA),strip.placement ="outside",strip.text =element_text(size =10))
2. Deterministic strength
Fitting models
Individual mountain models
Code
# iterate over each mountain in the dataset:# fit LM of each metric against linear and second-order polynomial terms of scaled elevation# extract model predictions for each mountainmetric_sess <-glue("SES.{metric_cols}")to_iter <-expand.grid(mountain = list_mount, metric = metric_sess)quad_reg_ses <-map2(to_iter$mountain, to_iter$metric, ~ { data_mount <- data_dryad %>%filter(Mountain == .x) model_formula <-as.formula(glue("{.y} ~ elevation_scaled + I(elevation_scaled^2)")) model_mount <-lm(model_formula, data_mount)# predict only for elevations present in each mountain to_pred <- data_mount %>%distinct(Mountain, elevation_scaled)# predict model_pred <-predict(model_mount, newdata = to_pred, type ="response") pred <-bind_cols(to_pred, pred = model_pred)tibble(mountain = .x, metric = .y) %>%left_join(pred, by =c("mountain"="Mountain"))}) %>%list_rbind() %>%relocate(metric, mountain) %>%# bring metric and mountain cols to the startmutate(metric =gsub("\\.", " ", metric)) %>%# remove period from metric namemutate(metric =factor(metric,levels =c("SES FRic", "SES FDis", "SES PD", "SES MPD")))# summarise original data to plot raw data pointsdata_dryad_ses <- data_dryad %>%select(Mountain, elevation_scaled, starts_with("SES")) %>%pivot_longer(cols =starts_with("SES"), names_to ="metric", values_to ="obs") %>%rename(mountain = Mountain) %>%mutate(metric =gsub("\\.", " ", metric)) %>%# remove period from metric namemutate(metric =factor(metric,levels =c("SES FRic", "SES FDis", "SES PD", "SES MPD")))
Global mixed-effects model
Code
quad_lmm_ses_output <-map(metric_sess, ~ { model_formula <-as.formula(glue("{.x} ~ elevation_scaled + I(elevation_scaled^2) + (1 | Mountain)")) model_ses <-lmer(model_formula, data_dryad, REML =FALSE) # because we use Likelihood Ratio Test# predict only for elevations present in each mountain to_pred <- data_dryad %>%distinct(elevation_scaled)# predict model_pred <-predict(model_ses, newdata = to_pred, type ="response", re.form =NA) pred <-bind_cols(to_pred, pred = model_pred) %>%mutate(metric = .x)# also generate model summary for Table 1 summ <-gen_mem_summ(model_ses) %>%mutate(metric = .x)# return bothlist(predicted = pred,summary = summ)})# model predictionsquad_lmm_ses <-map(quad_lmm_ses_output, "predicted") %>%list_rbind() %>%relocate(metric) %>%# bring metric col to the startmutate(metric =gsub("\\.", " ", metric)) %>%# remove period from metric namemutate(metric =factor(metric,levels =c("SES FRic", "SES FDis", "SES PD", "SES MPD")))# model summaryquad_lmm_ses_summ <-map(quad_lmm_ses_output, "summary") %>%list_rbind() %>%relocate(metric) %>%# refining formatting of labels and valuesmutate(metric =gsub("\\.", " ", metric), # remove period from metric nameterm =str_replace(term, "elevation_scaled", "Elevation")) %>%mutate(term =if_else(str_detect(term, "\\("),str_match(term, "\\(([^()]*)\\)")[, 2], term) %>%str_replace_all("\\^2", "<sup>2</sup>")) %>%mutate(term =str_replace(term, "Mountain", "Mountain (random effect)")) %>%pivot_wider(names_from ="metric", values_from ="value") %>%rename(Terms = term)
Lots of models!
# (elev. patt. for each mountain) + ((det. strength for each mountain) + (global det. strength))how_many_models <- (dim(to_iter)[1]) + (dim(to_iter)[1] +length(metric_sess))how_many_models
[1] 372
fig3_panel <-function(which_metric) { (# base plotggplot() +# raw data pointsgeom_point(data = data_dryad_ses %>%filter(metric == which_metric), col ="grey", alpha =0.5,mapping =aes(x = elevation_scaled, y = obs, group = mountain)) +# individual mountain model predictions (blue line)geom_line(data = quad_reg_ses %>%filter(metric == which_metric), col ="#00538F",mapping =aes(x = elevation_scaled, y = pred, group = mountain)) +geom_hline(yintercept =0) +# overall LMM prediction across all mountains (black line)geom_line(data = quad_lmm_ses %>%filter(metric == which_metric), linewidth =1.25,mapping =aes(x = elevation_scaled, y = pred)) +labs(x ="Elevation (scaled)", y = which_metric) +# match strip format in papertheme(axis.title.y =element_text(size =10)) ) %>%# include marginal histograms ggExtra::ggMarginal(type ="histogram", fill ="grey")}metric_sess_form <-gsub("\\.", " ", metric_sess)(wrap_elements(fig3_panel(metric_sess_form[1])) |wrap_elements(fig3_panel(metric_sess_form[2]))) / (wrap_elements(fig3_panel(metric_sess_form[3])) |wrap_elements(fig3_panel(metric_sess_form[4]))) &plot_annotation(tag_levels ="a", tag_prefix ="(", tag_suffix =")") &theme(plot.tag.position =c(0.05, 0.95)) # align tags with axis labels
Plotting everything
knitr::kable(quad_lmm_ses_summ)
Table 1: Parameter estimates (with SE) of mixed-effects models. The last row shows variances of the random effect term, Mountain.